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Abstract 

We propose a method to decompose the total energy of a supercell containing defects into 
contributions of individual atoms, using the energy density formalism within density functional 
theory. The spatial energy density is unique up to a gauge transformation, and we show that 
unique atomic energies can be calculated by integrating over Bader and charge-neutral volumes for 
each atom. Numerically, we implement the energy density method in the framework of the Vienna 
ab initio simulation package (vasp) for both norm-conserving and ultrasoft pseudopotentials and 
the projector augmented wave method, and use a weighted integration algorithm to integrate the 
volumes. The surface energies and point defect energies can be calculated by integrating the energy 
density over the surface region and the defect region, respectively. We compute energies for several 
surfaces and defects: the (110) surface energy of GaAs, the mono-vacancy formation energies of 
Si, the (100) surface energy of Au, and the interstitial formation energy of O in the hexagonal 
close-packed Ti crystal. The surface and defect energies calculated using our method agree with 
size-converged calculations of the difference between the total energies of the system with and 
without the defect. Moreover, the convergence of the defect energies with size can be found from 
a single calculation. 
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I. INTRODUCTION 



Total energy is one of the most important quantities in a solid state system, as it de- 
termines the stable configuration, and its derivatives provide other equilibrium properties. 
The assignment of energy to particular finite volumes can provide additional detailed infor- 
mation, such as defect formation energies. For example, the surface energy of a given facet 
of a crystal is meaningful for predicting the equilibrium crystal shape and preferred crystal 
growth directions, and should depend only upon the properties of the surface. The forma- 
tion energy of a point defect is important to understanding the phase stability, and should 
depend only upon the properties in the vicinity of the defect. However, density functional 
theory calculations^ provide only a single total energy for a given configuration, not a spatial 
partitioning of energy into additive contributions — instead, defect energies are determined 
as a difference of two of more separate total energy calculations. 

The energy density method^ can provide the formation energies for more than one point 
defect, surface or interface in a single calculation, as well as a picture of the distribution of 
the energy among the surrounding atoms. The energy density formula derived by Chetty 
and Martin is a reciprocal-space expression for norm- conserving pseudopotential 3 (NCPPs) 
with local density approximation^ (LDA), where the ion- ion interaction is treated in a 
similar manner as the Ewald sum.- The energy density is not unique since there are multiple 
definitions of the kinetic and Coulomb energy densities that integrate to the same well- 
defined total energy. The different definitions can be considered as gauge variability; defining 
a unique gauge-independent energy requires the identification of spatial volumes where the 
gauge differences integrate to zero. Chetty and Martin 2 ' 7 showed that the surface energy for a 
crystal can be calculated by an integral over a region to high-symmetry planes within the bulk 
of the crystal, where symmetry ensures that the gauge dependent terms integrate to zero. 
Therefore, two polar surface energies such as the (111) and the (111) surfaces of a zincblende 
semiconductor GaAs can be integrated independently in a single calculation. This compares 
with more complex approaches that use multiple wedge-geometries to extract polar surface 
energies for specific geometries.- Rapcewicz et a/.— followed the energy density formalism, but 
generalized the method to low-symmetry system such as the (0001) surface of GaN and the 
(0001) surface of SiC by introducing Voronoi polyhedra for the integration volumes; however, 
it should be noted that Voronoi polyhedra are gauge-independent integration volumes only 
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in specific situations. RamprasacU^ extended the application of energy density method from 
surfaces to point defects of metals and presented two applications on the monovacancy of 
Al and the (001) surface of Al. 

We reformulate the expression of the real-space energy density for the projector-augmented 
wave^i (PAW) method and the norm-conserving and ultrasoft pseudopotentials^ (USPP), 
grouping terms in a similar manner. We implement the energy density method in the Vienna 
ab initio simulation package^ 1 ^ (vasp). We also demonstrate the usefulness of assigning 
an energy to each atom using gauge-independent integration volumes. Excluding cases that 
can be determined by symmetry — such as the energy of an atom in the diamond structure 
that is half the total energy of a unit cell — the assignment is not unique. Nevertheless, there 
are two primary reasons for developing an approach based on gauge-independent volumes. 
One advantage is that it provides an automatic procedure to choose volumes that do not 
depend upon the choices of planes or polyhedra in the calculations mentioned above. The 
resulting surface and defect energies is the same as for the special cases above, but can also 
be used for general cases. In addition, the energy per atom partitions space to individual 
defect regions, and provides a measure of finite-size errors in a single calculation. Moreover, 
within an individual defect — e.g., a vacancy in silicon — the atomic energy shows the contri- 
bution of individual atomic relaxations to the total formation energy, which can be useful 
for understanding changes to the stability of individual defects and surfaces. 

The assignment of energy to an atom derives from Bader's "atoms in molecules" theory.— 
Although our work is different, one element is the same: the calculation of volume around 
each atom where the kinetic energy is unique. We extend this to find a charge neutral volume 
for a unique classical Coulomb energy. In Bader's work, a local form of virial theorem relates 
the local electron kinetic energy to the potential energy density, and defined atomic energies 
within each volume bounded by zero flux surface of the gradient of electron density. However, 
this results in charged units (so-called "Bader charges" ) which have long range forces. Thus 
it is arbitrary to assign the energy to one region when it is in fact shared by interacting 
regions. Also, the Bader approach does not consider the exchange-correlation energy of 
density functional theory. A recent paper— has reported a way to include such terms in 
a form of the virial theorem; they found different values from those given by the original 
method, making the applicability of a local virial theorem to density functional theory 
questionable. To overcome these difficulties, we instead define two integration volumes: one 
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based on kinetic energy, and another based on potential (electrostatic) energy. 

We derive the energy density methodology for PAW, and apply it to several defect types 
in solid state systems. Section [III reviews the energy density expression derived by Chetty 
and Martin for NCPPs, and derives the reformulated expression for the PAW method. The 
energy density contains a gauge-dependent kinetic energy density, a gauge-dependent long- 
ranged classical Coulomb energy density, a well-defined exchange-correlation energy density, 
and short-ranged terms grouped in an on-site energy for each ion. We also consider the 
method to integrate the gauge- dependent terms. Section [TTT] presents four applications: the 
(110) surface energy of GaAs, the mono- vacancy formation energy of Si, the convergence test 
of (100) surface energy of Au, and the interstitial formation energy of O in the hexagonal 
close-packed Ti crystal. These examples show the generality of the method, and the new 
information extracted such as convergence of formation energy with size and the spatial 
partitioning of defect energy. 

II. METHODOLOGY 

In density functional theory, the standard form for the total energy of a crystal is usually 
written in reciprocal space. There are many forms convenient for the total energy written in 
terms of the wavefunctions and/or eigenvalues.— For our purposes, we consider the expres- 
sions in real space where the Kohn-Sham wavefunctions are used only for the kinetic energy 
and for non-local terms in the pseudopotential. In the norm- conserving pseudopotential 
approximation, the total energy is, in atomic units (h = m e = e = 47r/eo = 1), 




nk 



+ 



£xc[p e (r)] 




(1) 
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where VVik(r) and / n k are the valence pseudo-wavefunction and the electron occupation 
number for the n th band, for wavevectors k within the first Brillouin zone, and with valence 
electron density p 6 (r). The first term is the independent electron kinetic energy. The second 
is the exchange-correlation energy E xc [p e (r)} = f dr p e (r)e xc (p e (r), | V(p e (r))|), where e xc 
is the exchange-correlation energy per electron; in the local density (LDA) or a general- 
ized gradient (GGA) approximation it is a function of the density or the density and its 
gradient. The fermion nature of many-body interacting electrons is approximated by this 
exchange- correlation potential. The third term is the energy due to the non-local part of 
the pseudopotential 



where is the Eih component of the non-local pseudopotential, with p£ the projection 
operator on angular momentum t. This term is site-localized (non-zero only within the core 
radius around a site) so that the total energy involves a sum over the sites p at position 
R M . The last three terms of Eqn. ([TJ) are the long-ranged Coulomb interactions. The fourth 
and fifth terms are the interaction of the electrons with the local ionic pseudopotential 
V^ oc (r — R M ) and with themselves that can be written as one-half the interaction with the 
Hartree potential V^(r) = / dr' ■ The sixth term is the valence charge of ion p at 
position R^ — the ion-ion Coulomb interaction energy that is the same as for point charges 
since the cores are assumed not to overlap. Finally, there can be non-linear core corrections 
not shown here, but which can be expressed in terms of E xc involving the core density 
similar to the PAW method. 

We use the PAW method^ for calculations, where the total energy has a similar form 
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+ E xc [p + p + p c ] 



(3) 



+ E K [p + p] 
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4=E^ii-~v%> 



+ £ xc [p 1 + Pc ] + £ H [p 1 ] 

+ / Vv\p Zc ](p X )dT 

' ~ 1 ~ (4) 



+ £ xc [p 1 + p + / 5 c ] + £ H [p 1 + p] 



for z,j = /me. Quantities with a tilde are obtained by pseudization, and a superscript 1 for 
quantities evaluated inside atom-centered spheres on a radial grid. For each atom-centered 
sphere, the pseudo-partial waves match all-electron partial waves \<pi) at the sphere 
boundary and outside the augmentation region. The smooth projector functions \pi) are 
dual to the pseudo-partial waves, and pjj = £) nk frikfynklPi) (Pjl^nk) are the occupancies of 
augmentation orbitals Then p is the soft pseudovalence electron density, p 1 and p 1 are 

the on-site charges (full and pseudized) localized around each atom, p is the compensation 
charge, p c and p c are the frozen core charges (full and pseudized), pz c and pzc are the sum 
of the nuclei pz and core charges (full and pseudized). The electrostatic interactions — 
electron-electron, electron-ion, and ion- ion interactions (last three terms in Eqn. (|3l)) — are 
collectively the "classical Coulomb" energy. The short-ranged terms for individual ions are 
-Eon- S ite = {E* — EV). We derive the total energy density for the PAW and pseudopotential 
methods as 



e(r) = t(r) + e cc (r) + e xc (r) + E on _ sitc 5(r - R M ) (5) 



and use gauge-independent integration over Bader and charge-neutral volumes to define 
atom-centered energies. 
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A. Kinetic energy density 



The kinetic energy density is gauge dependent and can be expressed as asymmetric or 
symmetric functional,— 

^ (a) (r) = -iE/^(r)V 2 ^ k (r) 

" k (6) 
^(r) = -^/n k |V^ k (r)| 2 . 

nk 

The difference between asymmetric and symmetric kinetic energy density is a gauge- 
dependent term proportional to the Laplacian of pseudo electron density, 

fM(r)-tM(r) = -^VV(r). (7) 

The integral of the two forms of kinetic energy density is equal when the gauge-dependent 
integral vanishes; e.g., for infinite or periodic systems. In Section Hi El we will integrate over 
a discrete set of atom-centered volumes where the gauge-dependent integrals also vanishes — 
hence, uniquely defined kinetic energies for atoms in a condensed system. Note that con- 
tinuous wavefunctions can have cusps in their gradient, thus the asymmetric form of the 
kinetic energy density can be ill-defined. Chetty and Martin chose the symmetric form for 
the kinetic energy density as it appears in the variational derivation of the Schrodinger's 
equation,™ and hence is a more fundamental quantity. However, the kinetic energy density 
is unique except for terms proportional to the Laplacian of pseudo electron density; if we 
integrate over volumes where the gauge-dependent term of Eqn. ([7j) (c.f. Section Hi El) , then 
either form of Eqn. ([6]) gives the same kinetic energy. For a planewave basis, the asymmetric 
kinetic energy density is well-defined everywhere — i.e., there are no cusps in the wavefunc- 
tion gradient — and is computationally less demanding to calculate. In the PAW method, 
the total kinetic energy density contains three terms 

t (a) ( r ) = 1(a) ( r ) + t l(a) ^ _ £l(a) ^ _ (g) 

The first term, (r) is also the first term in Eqn. ([3]) and is expressed by using the pseudo 
wavefunction. The on-site kinetic energies t 1(a) (r) and t 1(a) (r) are first terms in Eqn. (j3j), 
and are included in the short-ranged on-site energy E on _ sitc 5(r — R M ). 
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B. Classical Coulomb energy density 



The total classical Coulomb energy of a system with electrons and nuclei can be written 

as 



Z li Z v 
,,. „ /,> '"' 



(9) 



where p e (r) is the sum of soft pseudoelectron density p(r) and compensation charge p(r), fi 
and v are representing different nuclei and R^ v is the distance between two nuclei. There 
are various ways to calculate the electrostatic energy.~ In the Ewald method- terms are 
grouped with a Gaussian charge density around each atom so that the sum can be calculated 
by sums in both real and reciprocal space. However, this is not useful for constructing an 
energy density in real space. Instead, methods that involve only smooth densities for each 
ion can be used to construct expressions for the Coulomb energy that are expressed only in 
real space. 



1. Smeared ions 



We introduce a fictitious localized charge distribution p x ° c , which gives rise to a local 



pseudopotential V^ oc (c.f. Section F.3 of [17|) for ion \x as 



-—W c (r - R 



(10) 



The Coulomb interaction energy between two ions \x and v is 



- """" 



R 



I drp l ° c (r - R A1 )V r J oc (r - R v ) , 



(11) 



and the self energy on each ion is 



Ef = ^ / ^ P r(r)C C (r) ■ 



(12) 
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The total classical Coulomb energy of a system with electrons and nuclei can be written as 



2 J J |r-r'| 
+ J drp e (r) 5>^ c (r) 

+ V M. " (13) 



/,> '"' 



/ dr—\WV tot (v)\ 2 - V £ 

J 87T Z-/ - 



self 



with total classical Coulomb potential \^ tot (r) = V r H (r) + V loc (r). The Hartree, local, and ion- 
ion interaction terms (last three terms of Eqn. (j3J)) are combined into the classical Coulomb 
term. 

We use the Maxwell energy for the total electrostatic component of the energy density. 
As a fictitious charge density defined in Eqn. ffTU]) . the total neutral charge density p tot (r) = 
p e (r) + p loc (r). With this definition, the Maxwell form of the classical Coulomb energy 
density can be written as 

e Maxw e ii( r) = _L |w tot (r)| 2 . (14) 

Similar to the kinetic energy density, the classical Coulomb energy density is unique up to 
a gauge transformation. The asymmetric form of the classical Coulomb energy density is 

eg(r) = -iy"( r )VV tot (r) = V ot (r)p tot (r), (15) 

<57T I 

and the gauge-dependent term is the difference of Eqn. (Ill]) and Eqn. (115]) . 



J^IItA — ^Maxwell, 



1 



r) = -—V ■ [y tot (r)W tot (r)l . (16) 

As with the kinetic energy density, we can obtain a gauge-independent classical Coulomb 
energy as an integral over any volume bounded by a zero-flux surface of the gradient of the 
total Coulomb potential. 

The PAW method introduces the soft compensation-charge n, and the Hartree energy is 

E H =E H + (E H — E H ) 

=E H {p +p\+j2 W] - J2 2*1? + p] (17) 

The first term is related to the soft valence- charge density and the soft compensation-charge 
density, and is included in the classical Coulomb energy density. The last two short-ranged 



terms are related to the short-ranged on-site energy 

-^on-site "^7-0 ) *^ clfC the electron-ion 

interactions. 



2. Model charge density 




0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 



FIG. 1. Local pseudopotential and change density for PAW Ti (top), and compensating model 
potential and charge density (bottom). The PAW-GGA potential for Ti has a cutoff radius 1.22A; 
the charge density — given by the Laplacian of the potential — can have short- wavelength oscillations 
that are well-represented on a radial grid, but poorly represented on a regular Cartesian grid. To 
compensate for this, we add a smooth model potential and corresponding charge density to that 
matches outside of the cutoff radius. The model potential v model (r) is a smoothly varying long- 
ranged potential; the model charge density p model (r) is also a smoothly varying function for the 
radius r from to the cutoff radius r c , which integrates to the negative charge density of the local 
charge. 



In practice, the local charge density due to local pseudopotential can vary rapidly (c.f. Fig- 
ure HJ, which causes numerical errors in a real calculation. To improve our numerical accu- 
racy, we then introduce a model charge density p model ( r ) to solve this problem. The model 
charge density is chosen as obeying following constraints: a spherically symmetric functional, 
which is centered at each ion, zero beyond the cutoff radius of local pseudopotential, and 
normalized as negative to the local charge density within the cutoff radius. The total charge 
density can be rewritten as 

CIS) 



p tot (r) =p loc (r) + p model (r 



p-[T) - p (r) 
E W° c (r)+^ odel (r)]+^(r), 
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where pj^ c (r) + p™ odel (r) is a neutral and spherical charge density for each ion; 5p(r) is the 
difference between the valence electronic charge density and the model charge density. The 
asymmetric form of the classical Coulomb energy is 



E cc [p tot ] =E cc [p loc + p model ] 

+ J dr (V loc + V modcl )Sp 

+ dr5V s [5p] 5p 



(19) 



,-iself 



The first term is 

E cc [p loc + P m ° dcl ] = El° u c+modcl ( | | ) + 4 oc+model . (20) 

\1<V fl 

The electronic interaction between two neutral atoms is zero when there is no charge overlap, 
since all moments are zero for spherical charge distributions. Therefore, the first term in 
Eqn. ( 120]) is zero. Combining the second term with the self energy in Eqn. ffl2|) . we have 

£loc+modcl _ ^self = \ f dr (r) p^ (r) 

2J r (2i) 

+ J drVf c ( r )p™ Ac \v), 

which is a constant for each species of ions and can be canceled when studying the defect 
energies. Neglecting this constant term, the asymmetric form of the classical Coulomb 
energy density in Eqn. (fl9|) is 



-cc 



r) 



0° c (r) + V H (r) 

1 2 , (22) 

+ -\/ modcl (r)J [p e (r) - p modcl ( r ) 



where V loc (r), Vn(r) and p e (r) are already known in real space. Different model poten- 
tial v model (r) and model charge density p model (r) can be constructed as long as they obey 
above constraints. In this work, the model charge density is a polynomial functional with 
continuous zeroth-, first- and second-order derivatives at and r c ; for u = r/r c , 

, ^3 [1 - 10m 3 + 15m 4 - Qu 5 } :r<r c 

p model (r) = ^ 5*rg L (23) 

: r > r c 
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As shown in Figure [TJ the local charge density varies rapidly with respect to radius r, while 
the model charge density parametrized in Eqn. (1231) smoothly decays to zero as increasing r 
to r c . The corresponding potential is 



The model charge density gives faster numerical convergence on a regular spatial grid 
compared to the rapidly varying local charge density. This model charge density has been 
tested by calculating the surface energy based on a Si bulk (8 atom) and a Si (111) slab (16 
atom). The energy difference between the values calculated from total energy calculation 
and from energy density integration is about 7meV. It has also been tested on O atoms, and 
O2 molecules for various grid sizes, where we had a convergence problem by using a rapidly- 
varying local charge density. For a wide range of grid sizes, the total energy calculations 
converge to a precision of lmeV, while the difference between results calculated from those 
two methods is up to 0.4eV. However, the energy difference can be reduced below lmeV 
with the smooth model charge. 

C. Exchange- correlation energy density 

The exchange-correlation energy of many-body interacting electrons can be expressed 
in terms of an exchange-correlation hole which tends to be localized around each electron. 
In density functional theory it is usually treated as a function of the local density and its 
gradients, which is determined by the choice of exchange-correlation functional. For both 
the local density approximation (LDA) and generalized gradient approximation (GGA), the 
gauge-independent exchange-correlation energy density is 



where e xc is exchange- correlation energy per electron. 
D. On-site energy 

The last term of the energy density in Eqn. fl5]) is short-ranged. For the PAW method, 
the on-site energy for each ion is composed of kinetic energy, exchange-correlation energy, 




(24) 




(25) 
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and Coulomb energy including electron-electron and electron-ion interactions in the aug- 
mentation region, and is 

E on . site = (El - £>(r - R„) , (26) 

with the on-site energies El and E^ expressed in Eqn. (jlj). In practice, we calculate 
and E^ for each ion using radial grids. For NCPPs and USPPs, this short-ranged term 
corresponds to the non-local pseudopotential energy. For USPPs, 

E t = E/ dr ^) ( E D i?w<&i ) ^k(r) , (27) 

nk \ ij / 

where the coefficients D^ n and the projector functions vary depending on the atomic 
species. For NCPPs, 

< = EE / <Mk(r)V#(|r-R M |)p £ ^(r). (28) 
E. Gauge dependence and uniqueness 

We have two energy density terms to be integrated which are gauge dependent: the kinetic 
energy density and the classic Coulomb energy density. Defining a gauge independent energy 
requires integrating these energy densities over volumes to cancel out any gauge dependence. 
Previously, Chetty and Martin 2 - integrated over Wigner-Seitz cells in a supercell; Rapcewicz 
et ai— constructed a Voronoi polyhedron for each comprised atom. However, these volumes 
are not the general solution to removing gauge dependence. For the kinetic energy density, 
the gauge dependence, Eqn. ([7]), is proportional to the Laplacian of electronic charge density. 
Hence, we integrate over a volume where the gradients of electron density has zero component 
along surface normal direction n, Vp(r) ■ n = — the zero-flux "Bader" volume. 15 For the 
classical Coulomb energy density, the gauge dependence, Eqn. f JTB"]) . is proportional to the 
Laplacian of the potential. Hence, we integrate over a volume where the electrostatic field 
has zero component along surface normal direction n, VV"(r) • h = — the zero-flux charge- 
neutral volume. 

In this work, we construct two different volumes: the Bader volume is used to integrate 
kinetic energy density and exchange-correlation energy density, and the charge-neutral vol- 
ume is used to integrate classical Coulomb energy density. Each of these volumes is "atom- 
centered" — that is, it contains one atom somewhere in the volume — and they each partition 
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space: the union of all volumes is the total supercell volume, and the intersection of any 
two volumes is zero. We define these volumes on the same regular spatial grid used to rep- 
resent the charge density and energy density terms. Accurate definition of the volumes and 
integration uses a weighted integration scheme^ that has quadratic convergence in the grid 
density. Finally, we use the integral over the gauge- dependent kinetic (Eqn. ([7])) and classi- 
cal Coulomb (Eqn. (TTBjl ) terms as an estimate of the integration error for atomic energies. 
This error estimate has a sign, so it is possible for the magnitude of the error of the energies 
of neighboring volumes A and B to be greater than the magnitude of error integrating over 
AuB. 



F. Summary 

We summarize the procedure of calculating atomic energy using energy density method 
in Table [B The total energy density of Eqn. §5§ contains a gauge-dependent kinetic energy 
density, a gauge- dependent classical Coulomb energy density, a gauge-independent exchange- 
correlation energy density, and a short-ranged non-local energy. The well-defined atomic 
energy can be calculated with two different integration volumes. The Bader volume is 
employed for the integral of the kinetic energy and the exchange-correlation energy, and the 
charge neutral volume for the integral of the classical Coulomb energy. 



III. APPLICATIONS 



To verify our implementation of the energy density method and highlight the new in- 
formation it reveals, we perform DFT calculations with vasp on the GaAs(llO) surface, Si 
monovacancy, Au(100) surface and O interstitial in Ti. We integrate the energy densities 
around the defect regions, and compare the integrated defect energies with values given 
by total energy calculations and experiments. Finally, the convergence of the atomic en- 
ergy to bulk values within a single calculation shows the convergence (or lack of) for each 
calculation. 
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e(r) = f(r) + e cc (r) + e xc (r) + S on _ Bito <5(r - R M ). g) 

1. Kinetic energy density 

* (s) (r) = |E„ k /nk|V^ nk (r)| 2 . 

t(»)( r ) = —I £ nk / nk ^ k (r)V 2 ^ nk (r). © 
i(">(r)-fW(r)=-ivV(r). 
Construct zero-flux volume O p where Vp e (r) ■ h = 
The bounded volume integral, T = J n £(r). 

2. Classical Coulomb energy density 

e Maxwoll( r ) _ ^_| W tot( r )|2_ (T3J 

e cc( r ) = |V tot (r)p tot (r). (13 

= [V loc (r) + |V H (r) + |y model (r)][p e (r) - p modol (r)]. 
e^(r) - e^9 xwell (r) = _^_v ■ (V tot (r)VV tot (r)). g6} 
Construct zero-flux volume f2y where VV tot (r) ■ n = 
The bounded volume integral, i?cc = / n ecc(r). 

3. Exchange-correlation energy density 

exc(r) =p e (r)exc(p e (r)). HD> 
The bounded volume integral, -Exc = Jq e x c(r). 

4. On-site energies 

PAW: Somite = - El). lH 

USPPs: S»> = £ nk / rfr^* k (r)(Ey f {flftX^ l)ik(r). 123 

NCPPs: = E nk E* / rfr</> k (r)V$(|r - R M |)p^„ k (p). 121 
The atomic energy: E = T + E C c + E xc + Son-site- 

TABLE I. Summary of the energy density formulae for PAW, USPPs, and NCPPs methods and 
the procedure to calculate atomic energy using the energy density method. 

A. GaAs(llO) surface 

The GaAs(llO) surface contains equal numbers of Ga and As atoms: a stoichiometric or 
non-polar surface. The surface energy 7 sur f of a stoichiometric slab is 

7surf = ^ ( -Eslab — N s \ ah jy bUlk j ) (29) 



bulk 



15 



for surface area A, where -E s i a b is the total energy of a GaAs slab with iV s i a b pairs of GaAs 
atoms, and -Etmik is the total energy of GaAs bulk with iVtmik pairs of atoms. Our DFT 
calculations are performed with the PAW method,— with the local density approxima- 
tion (LDA) 4 -^ for the exchange-correlation energy The valence configurations for Ga is 
[Ar]3d w As 2 Ap 1 with cutoff radius l.OlA, and As is ([Ar]3c/ 10 )4s 2 %> 3 with cutoff radius l.llA; 
this requires a plane-wave basis set with cutoff energy of 650eV. This gives a lattice constant 
of 5.6138A for zincblende GaAs, compared with the experimental lattice constant of 5.65A. 
The supercell contains 11 layers of atoms with a pair of GaAs atoms on each layer, and a vac- 
uum gap of 8A to prevent the interaction between slabs under periodic boundary conditions. 
We use Monkhorst-Pack k-point meshes 20 of 8 x 8 x 8 for bulk eight-atom cells, and 8x8x1 
for the slab supercell; Brillouin-zone integration uses Gaussian smearing with k^T = O.leV 
for electronic occupancies, and the total energy extrapolated to k^T = OeV. We represent 
the charge density and compute energy densities on a grid of 84 x 120 x 560. Geometry is 
optimized to reduce forces below 5meV/A. This gives a surface energy of 50meV/A 2 ; this 
agrees with Moll et aVs valued of 52meV/A 2 , Qian et aUs value 2 ^ of 57meV/A 2 , Choud- 
hury et a/.'s 2 ^ LDA value of 50meV/A 2 , and the experimental value 2 - 4 of 54 ± 9meV/A 2 . 

Figure |2] shows the energy change from bulk for each layer by integrating over volumes 
that eliminate gauge dependence. The change in energy shows differences from bulk that 
are mainly confined to the first two layers; the bulk-like response of the interior layers — not 
just for the total energy, but also the individual contributions to the energy. Determining 
the size- convergence of a surface calculation with total energy alone requires computing 
surface energies for multiple sizes; in our case, the bulk-like behavior of our center layers 
indicates a small finite-size error without requiring multiple size calculations. We can inte- 
grate the surface energy by adding the energies from the first two layers; our surface energy 
is 51 ± lmeV/A 2 , which agrees well with the total-energy calculation of surface energy. The 
error estimate is specifically for the integration error over the Bader and charge-neutral vol- 
umes. Note also that we can compute the energy of each surface independently; for surfaces 
with different chemistry, this allows for two surface energies to be calculated from a single 
supercell. 

Figure |3] shows the Bader and charge-neutral volumes for Ga and As atoms in a (110) 
plane of GaAs. The Ga and As atoms all lay in a plane, and the intersection of the surfaces 
show the difference between the two atom-centered volumes. The Bader volumes have zero 
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FIG. 2. Atomic energy distribution on GaAs(llO) slab. The supercell contains 11 layers of GaAs. 
The energy density integrated over each atomic layer divided by surface area gives the energy per 
layer referenced to bulk value A£j aycr . The atomic integration errors are smaller than lmeV/A 2 . 
Surface energy is confined to first two layers. The bulk-like behavior of the center layers indicates 
the sufficient thickness of slab calculation. The individual energy term contributions to each layer 
is shown in bottom plot. All the energy terms are bulk-like for the center layers of the slab, not 
just the sum. At the surface, kinetic energy decreases as the valence charge density spreads out 
into the vacuum. 
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FIG. 3. (Top) Bader volumes and (bottom) charge-neutral volumes for Ga (purple) and As (yellow) 
in a (110) plane of GaAs. These integration volumes define the unique atomic kinetic and exchange- 
correlation energies, and classical Coulomb energies. Each volume contains a single atom, but the 
two volumes are different for the same atom. 

flux of the gradients of charge density through their surfaces, and are used to integrate a 
unique kinetic and exchange-correlation energy. The charge-neutral volumes have zero flux 
with the total electrostatic field, and are used to integrate a unique classical Coulomb energy. 
These volumes are different also from the Voronoi volumes around each atom. The atomic 
volumes, like the individual components of energy, become bulk-like in the center of slab. 
Atoms at the free surfaces have volumes that extend into the vacuum. Besides the different 
surfaces, the Bader volumes of As are larger than the As charge-neutral volumes. 



B. Si monovacancy 

The monovacancy in bulk Si is a simple point defect in a semiconductor, which has been 
studied theoretically and experimentally. From total energies, the formation energy of a 
vacancy AH V is 

AH V = - ^E N , (30) 

where E^ -1 and E N are the total energy of the N — 1 and A" atom supercells with and 
without one vacancy. Wright^ performed LDA- and GGA-PBE^ calculations in 215, 51 1-, 
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and 999- atom supercells to get formation energies of 3.53eV, 3.49eV, and 3.47eV with LDA 
and 3.66eV, 3.63eV, and 3.62eV with GGA. Puska et al— performed LDA calculations in 
31-, 63-, 127-, and 215-atom supercells to get formation energies of 3.98eV, 3.42eV, 3.44eV, 
and 3.31eV. Experiments have found a formation energy of 3.6 ± 0.2eV.— 

Our DFT calculations are performed with the PAW method with the generalized gradient 
approximation (GGA) of Perdew and Wang (PW91)^ for the exchange- correlation energy. 
The valence configurations for Si is [Ne]3s 2 3p 2 with cutoff radius l.OlA; this requires a 
plane-wave basis set with cutoff energy of 417eV. This gives a lattice constant of 5.4674A 
for diamond Si, compared with the experimental lattice constant of 5.43A. The 2x2x2 
simple cubic supercell with a vacancy contains 63 atoms. We use a 4 x 4 x 4 Monkhorst-Pack 
k-point mesh; Brillouin-zone integration uses Gaussian smearing with k-^T = 0.15eV, and 
the total energy extrapolated to k^T = OeV. We represent the charge density and compute 
energy densities on a grid of 200 x 200 x 200. Geometry is optimized to reduce forces below 
5meV/A. This gives a formation energy of 3.65eV. 

Figure H] shows the energy change from bulk for shells surrounding a Si vacancy. The 
primary contribution to the vacancy formation energy comes from the first five shells, be- 
coming bulk-like at larger distances. Summing the atomic energies up to the fifth shell gives 
a formation energy of 3.57 ± 0.05eV, which is similar to the total energy calculation. The 
importance of the fifth shell over the third and fourth shells can also be seen in charge 
disturbances from a vacancy. Kane 30 showed charge disturbances around a Si monovacancy 
out to the 27 th shell, while the first two shells contribute 60% of the charge disturbance. The 
most striking future was charge concentrates on the {110} planar zigzag chains of atoms, 
such as [000], [111], [220], [331], [440], as so on. Twelve such chains exist by symmetry. 
After reaching the fifth shell at |(331), the charge decays monotonically along the coplanar 
chains. Kane connected this result to the importance of the fifth-neighbor interaction in the 
valence force model 3 ^ of covalent phonon spectra. The valence force model is an empirical 
model connecting force constants to the electronic configuration. The fifth-neighbor interac- 
tion is proportional to r 2 A(pA(p' from changes in bond angles <p and ip' along a zigzag chain. 
The fifth neighbor has a stronger interaction with the bond-bending than third and fourth 
neighbors; we see a similar change in energy for the vacancy. 
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FIG. 4. Si monovacancy (red) in 2 x 2 x 2 simple cubic supercell with atomic energies. The first 
five shells are § (111) orange, § (220) yellow, § (113) green, f (004) blue, and f (331) violet. Only 
three shells have displacements greater than O.OlA: first shell with 4 atoms relax inwar d by 0.17A; 
the second shell with 12 atoms, relax inward by 0.05A; and the fifth shell relax inward by 0.02A. 
These three shells form a zigzag chain (red) from the vacancy with a strong interaction. The 
atomic energy confirms this interaction with energies of 560 ± 5meV (first), 49 ± 12meV (second), 
47 ± lmeV (fifth); compared with 16 ± 7meV (third) and — 4 ± 3meV (fourth). As with a free 
surface, the kinetic energy drops close to the vacancy due to decreasing of valence charge density. 
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C. Au(lOO) surface 



Au(100) is a low- index metallic surface without a reconstruction. To calculate the surface 
energy with Eqn. (|29[) . the thickness of the slab must be increased until the surface energy 
7surf converges. Our DFT calculations are performed with the PAW method with the gen- 
eralized gradient approximation (GGA) of Perdew, Burke, and Ernzerhof (PBE) 26 for the 
exchange-correlation energy. The valence configurations for Au is [Xe]6s 1 5<i 10 with cutoff 
radius 1.32A; this requires a plane- wave basis set with cutoff energy of 400eV. This gives a 
lattice constant of 4.17lA for FCC Au, compared with the experimental lattice constant of 
4.08A. Supercells range from 4- through 7-layers of atoms with one Au atom on each layer, 
and a vacuum gap of 10. 5A to prevent the interaction between slabs under periodic boundary 
conditions. We use Monkhorst-Pack k-point meshes of 13 x 13 x 13 for bulk four-atom cells, 
and 13 x 13 x 1 for slab supercells; Brillouin-zone integration uses the Methfessel-Paxton 
method^ 2 - with k-^T = 0.2eV for electronic occupancies, and the total energy extrapolated 
to k B T = OeV. We represent the charge density and compute energy densities on a grid 
increasing from 60 x 60 x 350 to 60 x 60 x 460 for 4- to 7-layer supercells. Geometry is 
optimized to reduce forces below 5meV/A. This gives surface energies of 53.3, 52.7, 52.5 to 
52.3meV/A 2 ; this agrees with Zolyomi et a/.'s valued of 54meV/A 2 . 

Figure |5] shows the energy change from bulk value for each layer of atoms. Energy density 
integration of left (or right) 2 layers of 4-layer slab gives a surface energy of 52 ± lmeV/A 2 ; 
5-, 6-, and 7-layer slabs all give a surface energy of 50 ± lmeV/A 2 . Although the calculated 
surface energy of 4-layer surface energy is close to the 5-, 6- and 7-layer 's values, the atomic 
energy distribution shows the center layers of 4-layer slab have not reached the bulk-like 
behavior. Unlike the convergence test of the traditional total energy calculation, the required 
thickness of a slab can be directly determined by observing the variation of atomic energy 
from the bulk value. 

D. HCP Ti with O interstitial 

Finally, we consider the formation energy of an oxygen interstitial in the octahedral site 
of HCP titanium. The formation energy is 

Ei'' = E(Ti + Oi) - E(Ti) - l -E{0 2 ), (31) 
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FIG. 5. Atomic energy distribution on Au(100) slab. The atomic integration errors are smaller 
than lmeV/A 2 . The bulk-like behavior of the center layer(s) of 5-, 6-, and 7-layer slabs indicates 
the sufficient thickness of slab calculation. 



where E(Ti + O;) and E(Ti) are the total energy of relaxed supercells with and without 
an oxygen atom, and E(C>2) is the total energy of oxygen molecule. Our DFT calculations 
are performed with the PAW method with the GGA-PW91 for the exchange-correlation 
energy. The valence configurations for Ti is [Ne]3s 2 3p 6 4s 2 3<i 2 with cutoff radius 1.22A, and 
O is [He]2s 2 2p 4 with cutoff radius 0.80A; this requires a plane-wave basis set with cutoff 
energy of 500eV. This gives a lattice constant of a =2.933A, c =4.638A, and c/a =1.581 
for HCP Ti, compared with the experimental lattice constant of a=2.95lA, c =4.684A, 
and c/a =1.587.~ The supercell contains 96 Ti atoms (4x4x3) and 1 O atom. We use a 
2x2x2 Monkhorst-Pack k-point mesh; Brillouin-zone integration uses the Methfessel-Paxton 
method with k^T = O.leV for electronic occupancies, and the total energy extrapolated to 
ksT = OeV. We represent the charge density and compute energy densities on a grid of 180 x 
180 x 216. Geometry is optimized to reduce forces below 20meV/A. This gives an oxygen 
interstitial formation energy of -6.19eV, with a nearest-neighbor distance between Ti and 
O of 2.08A. Hennig et al— performed GGA-PW91 calculations using ultrasoft Vanderbilt- 
type^ pseudopotentials in the same supercell to get a formation energy of -6.12eV and 
nearest-neighbor distances of 2.06-2.09A. 

Figure [6] shows the calculated Ti atomic energy change from bulk value for each shell. 
The change in energy shows the differences from bulk are mainly confined to the first two 
shells, with bulk-like behavior for shells further away from the oxygen atom. The energy 
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FIG. 6. HCP Ti 4 x 4 x 3 supercell with an O interstitial (red) in an octahedral site. The formation 
energy is confined to the first two Ti shells (yellow); away from the oxygen atom, the charge density 
and energy density of Ti atoms experience Friedel-like oscillations. There is a charge transfer of 
1.45e to the interstitial O atom. 



density and charge density oscillates and decays away from the interstitial. They peak 
at the 6th shell and the 13th shell with a wavelength of 1.9A. Weiss's Compton profile^ 
measured the Fermi momentum of Ti as 1.08 ± 0.06a.u., which corresponds to a Friedel 
oscillation wavelength of 1.5 A. Jepson's^ earlier calculation using linear muffm-tin-orbital 
method obtained the Fermi energy of Ti as 0.667Ryd which corresponds to the Friedel 
oscillation wavelength of 2.0A. Adding the atomic energy change of first two shells from the 
O interstitial and the atomic energy change of O atom, we obtain the interstitial formation 
energy of —6.13 ± O.OleV, which agrees with the total-energy calculation. 
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IV. CONCLUSIONS 



We implement the energy density method for PAW and USPPs for the planewave DFT 
code VASP; and analyze surface energies from the energy density in the surface region; 
and vacancy and interstitial formation energies from the energy density in the point defect 
region. The method can be applied to surfaces and defects in a variety of system, and 
produces defect formation energies comparable to well-converged total energy calculations. 
Furthermore, the energy density determines the distribution of energy near the defect or 
surface, and the sufficiency of a supercell without a separate convergence test. It can also 
give separate defect formation energies from a single supercell calculation. 
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